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NUMMARY 


Numerical solutions of the interacting laminar boundary- 
layer equations are presented for two symmetric airfoils at zero 
incidence: the NACA 0012 and the NACA 66^-018 airfoils. The 
potential flow was computed using Carlson’s code, and viscous 
interaction was treated following a Hilbert-integral scheme due 
to Veldman. Effects of various grid parameters are studied, 
and pressure and skin-friction distributions are compared at 
several Reynolds numbers. For the NACA 0012 airfoil, Reynolds 
number is varied from a value just below separation (R N = 3000) 
to a value for which extensive separation occurs (R^ = 100,000) . 
For the 66^-018 airfoil, results are given at intermediate 
values (R = 10,000 and 40,000). The method fails to converge 
for greater values of Reynolds number, corresponding to the 
development of very thin well-separated shear layers where 
transition to turbulence would occur naturally . 
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INTRODUCTION 


A major problem in computational aerodynamics is the predic- 
tion of maximum lift coefficient for two-dimensional wing sections. 
The numerical solution of this problem would involve computation 
of the flow past airfoils at large angles of attack, including 
conditions of separation and stall. At large Reynolds number the 
boundary— layer model is appropriate, provided viscous interaction 
is included; in fact, separation bubbles of moderate extent have 
been analyzed on the basis of such a model. The long-term goal 
of this study is to carry out such computations for airfoils at 
moderate angles of attack, with the hope that a computation of this 
type would provide aerodynamic data up to inception of massive 
separation. However, the present study has been restricted to 
zero angle of attack, with emphasis on the complications of laminar 
separation arid wake flow that arise with that condition. 

The computer program developed for this study is based on the 
laminar-flow condition. This was desirable for two reasons: 

(1) leading-edge separation, which can limit the maximum lift, is 
usually of laminar type, and (2) the basic interaction phenomenon 
is the same for laminar and turbulent flow, while the equations 
of motion are simpler in the former case. In addition, laminar 
flow solutions are of interest for their own sake . For higher 
Reynolds numbers, an appropriate turbulence model can be 
incorporated . 

Various computer codes are available for the potential— flow 
portion of the computations. A code developed by Carlson (Ref. 1) 
for airfoil analysis and design was selected for adaptation to our 
use; it is well-documented and structured in a way that readily 
permits interaction with a boundary-layer code in both direct and 
inverse modes of solution. Since our current interest is in low- 
speed flow, a simplified version of Carlson's code was programmed 
for the special case of incompressible flow. In addition, the 
multigrid technique was incorporated to accelerate convergence of 
the potential-flow computations . Our implementation of the 
multigrid procedure, as applied to Carlson's program, is described 
in Appendix B . 
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A new boundary-layer code was developed based on a stream- 
function formulation that permits a natural treatment of inter- 
action. Direct, inverse, and interaction modes all are included 
in the solution procedure, with switching between modes permitted 
at any streamwise station. Details of the method are given in the 
body of this report. 

An attempt was made to couple the potential-flow and bourldary- 
layer codes and operate each in the inverse mode. However, a 
converged solution was not achieved by this method. Instead, an 
alternative fully interacting scheme was developed successfully . 
This method and results obtained with it are discussed in the 
following sections . 


INTERACTION THEORY 


In the classical boundary-layer theory proposed by Prandtl, 
the pressure gradient acting on the boundary layer is pre- 
scribed by the inviscid (potential flow) solution and the 
boundary-layer properties, such as the displacement thickness, 
result from the solution. Prandtl (Ref. 2) also proposed that 
the pressure distribution could be determined to higher-order 
accuracy by recalculating the potential flow, accounting for 
the displacement thickness of the boundary layer. This scheme of 
accounting for the viscous effect on the potential flow is 
known as weak-interaction theory. Unfortunately it fails when 
boundary-layer separation occurs, since Goldstein (Ref. 3) has 
shown that the boundary-layer solution at separation is singular 
when the pressure gradient is prescribed, so that an accurate 
computation must break down at the separation point. Catherall 
and Mangier (Ref. 4) showed that this difficulty can be overcome 
by treating the boundary-layer problem as an inverse problem, thus 
solving for the pressure distribution which would produce a 
specified displacement thickness. Of course, in practice neither 
pressure gradient nor displacement thickness are known a priori. 

The successful treatment of flow separation as an inverse 
problem, contrasted with the failure of the classical theory, 
indicates that separated flows are of the strong-interaction 
type; i.e., the equations of the external potential flow and of 
the boundary layer must be coupled and solved simultaneously. 

A qualitative theory of viscous interaction based on this concept 
was presented by Crocco and Lees (Ref. 5). The idea also appears 
in the early work of Lighthill (Ref. 6), and a rational theory has 
evolved in more recent times in the independent works of Neiland 
(Ref. 7)» Stewartson and Williams (Ref. 8), and Messiter (Ref. 9). 
This theory, termed the triple-deck by Stewartson, is valid as an 
asymptotic solution of the laminar Navier-Stokes equations in the 
limit of infinite Reynolds number. (We note that Prandtl's 
boundary-layer theory is a valid limit solution for unseparated 
flows in the same sense.) 
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The first complete problem to be solved in triple-deck 
theory was that of the viscous interaction at the trailing edge 
of a flat plate*. This problem is an example of one for which 
the classical weak-interaction theory fails. The interaction 
is caused by the discontinuous change from a no-slip condition 
on the plate to a no-stress condition on the wake centerline. 

The resulting solution of the classical boundary-layer theory 
exhibits a singularity in the slope of the displacement surface 
for x = 0 (x = 0 at the trailing edge), and a corresponding 
singularity then is produced in the lnviscid pressure distribu- 
tion if the computation is pursued in the classical manner. In 
triple-deck theory, the viscous lower-deck solution is determined 
simultaneously with the inviscid upper deck and no singularity 
develops. The solution was obtained first by Jobe and Burggraf 
(Ref. 10), and confirmed by Veldman and Van de Vooren (Ref. 11) 
and by Chow and Melnik (Ref. 12). A summary of the results is 
given by Figure 1. The skin friction shown there is normalized 
by the Blasius (non-interacting) value; the effect of interaction 
is seen to raise the value by about one third of the Blasius 
value at the trailing edge. Correspondingly, the pressure on the 
plate is reduced below the freestream value, and then quickly 
rises in the wake to a level above the freestream value with a 
slow decay downstream. 

Chow and Melnik (Ref. 12) went on to consider the flat plate 
at angle of attack. The triple deck then substitutes for the 
Kutta condition in classical airfoil theory, and in fact the so- 
lution gives the viscous correction to the ideal lift of the air- 
foil. For fixed (but large) Reynolds number, viscous loss of 
lift is a relatively weak function of angle of attack for small 
a, but appears to increase catastrophically as a approaches 
the stall limit. Actually the computation scheme failed before 
the stall limit was reached, but an estimate was made by extrapo- 


"Stewartson and Williams (Ref. 9) solved for the upstream eigen- 
functions in supersonic flow, but did not match these to particu- 
lar downstream disturbances. 


6 


lation of the computed results to the point at which separation 

first occurred at the trailing edge. Thus the first rational 

prediction of C. (maximum lift coefficient) was made on the 

niax 

basis of triple-deck theory. Of course for the flat plate in 

laminar flow, C. is likely to be determined by leading-edge 

max 

separation, whereas the analysis of Chow and Melnik is restricted 
to the vicinity of the trailing edge. Nevertheless, a major step 
in the theory has now been taken. 

In reality, the Reynolds numbers for laminar flow are not 
high enough for triple-deck theory to be quantitatively accurate. 
Consequently for practical calculations the interacting boundary- 
layer theory is more appropriate. This theory was first thought 
to be merely an approximate model of the flow processes, but it 
is now known that it is exact in the limit of infinite Reynolds 
number in the same sense as is triple-deck theory (Burggraf, et al. 
Ref. 26). The main source of error in the triple-deck appears 
to be the treatment of the upstream boundary conditions. Because 
of its asymptotic nature, the longitudinal length scale approaches 
zero, and hence, the upstream boundary condition (matching to the 
boundary layer) is applied infinitely far upstream on the triple- 
deck scale. In contrast, in the interacting boundary- layer model, 
the upstream condition is that the boundary layer originate at the 
leading edge. The simplification in the triple-deck case is that 
Reynolds number can be scaled out of the problem, whereas it must 
be specified as a parameter in interacting boundary-layer theory. 
This complication is more than made up for by the improvement 
in accuracy for practical applications. Experience with the 
interacting boundary-layer theory indicates that accurate re- 
sults can be obtained for flows that include regions of separ- 
ation having moderate transverse extent. The model is likely 
to fail when breakaway occurs, with its resultant catastrophic 
reduction of lift. Thus even if breakdown of the theory occurs, 
it would be an indication of maximum lift conditions. 
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The computations described above were based either on a time- 
marching or similar procedure (Rizzetta et al, Ref. 13) or an in- 
verse method of solution (Jobe and Burggraf, Ref. 10; Veldman and 
Van de Vooren, Ref. 11; Chow and Melnik, Ref. 12). V/ith the latter 
method, the potential-flow problem is solved for the shape of the 
displacement surface for specified pressure distribution, and con- 
versely for the boundary-layer problem.. (The direct method of 
solution is unstable.) Both of these methods are slow, requiring 
severe under-relaxation or small time step to converge; thus long 
computer runs are required to achieve the steady-state solution. 
Recently Veldman (Ref. 1 ; 1) has presented an alternative method, 
which appears to be significantly faster than either of the above 
methods. In his method, the potential flow and boundary layer are 
coupled on a given vertical line and solved together, station by 
station. Upstream influence is accomplished by iteration on the 
whole flowfield solution. The method was shown to work well for 
both unseparated wake flows and for separated flow over contoured 
surfaces. Consequently, it would appear to be a good method for 
separated flow past finite airfoils. As will be seen, however, 
the efficiency of the method deteriorates when the combination of 
flow separation and wakes occurs . 
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FORMULATION 


Viscous-Interaction Condition 

The inverse method Invokes viscous interaction by a back-and- 
forth iterative transfer of information between the potential flow 
and boundary-layer codes. In contrast, Veldman*s (Ref. 14) idea 
was to compute the boundary layer and potential flow simultane- 
ously station-by-station so that the two flows fully interact at 
all stages of the iteration process. In Veldman's work, the 
potential flow was represented by the Hilbert integral of line- 
arized theory, relating pressure to the slope of the displacement 
surface. This formulation is strictly valid only for small per- 
turbations of the displacement surface from a true plane. However, 
the Idea is readily adapted to a potential flow computation pro- 
cedure like that of Carlson. 

According to linearized theory, the interaction law can be 
expressed In terms of a Hilbert integral as 

oo 

, l| v e (x) dx 


J 

_00 


where (u e ,v e ) are the (x,y) components of velocity at the edge of 
the boundary layer. In terms of the surface ordinate y g and the 
displacement thickness of the boundary layer 

v e (x) = u » 55 <y s + S,) 


The inviscid surface speed u g results from the above Hilbert 
Integral if 6* is set equal to zero. Hence the deviation of u e 
from u is given by „ 

s r' 


u 


e 


U c 

u + — 

S 71 


J 


d6* dx 
d3T x-X 


( 1 ) 
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Thus In (1) u can be calculated by any potential-flow method, 
s 

while the displacement correction due to the boundary layer is 
the given by the Hilbert Integral of linearized theory. In the 
present study, u has been evaluated using the Carlson code. In 
the discretized version the Hilbert integral is truncated at 
finite limits, resulting in a Cauchy integral; then, using the 
trapezoidal Integration rule, the integral Is replaced by a 
finite sum. Denoting values at the i-th grid point along the 
x-axis by the subscript i, we have 

Ue i = Us i + ° lj6j (2) 


Veldman gives the values of the c^ corresponding to a uniform 
grid; we have generalized his derivation to allow for a non-uniform 
grid, as presented in Appendix A. 

In Eq. (2) neither u nor 6* are known a priori, but are 
related by the solution of the boundary-layer equations, which 

are solved sequentially at successive stations. Thus at a given 
station x 1 , the values of 6*, for j < i are known (at least In the 
current approximation) while those for J > i are known only In the 
preceding approximation. For the values at j = i, Eq. (2) Is 
used In the form 


u 


u 


'ii 


6 > 



M 

y~ c, . 6 # 
~R iJ j 

j/i 


( 3 ) 


together with the no-slip condition at the surface 

u = v = 0aty = y s 

The boundary-layer equations and the above interaction condition 
are solved iteratively by marching from x q to x^ in each iterate 
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with the unknown terms for j M on the right side of Eq. (2) 
evaluated from the known solution for the previous iterate. In 
our application of Eq. (3), a relaxation factor r v was introduced 
to allow more flexibility with the calculation; letting the super- 
script denote the iteration number, we have 




* (n) 




(n-1) 



6 


*(n-l) 

i 




Thus r v ^0 produces the direct problem, with Goldstein singularity 
at separation, r v = 1 corresponds to Veldman's approach, and r v -*■ 00 

yields the inverse problem. Usually the value r y = 1 is satis- 
factory, but for the higher values of Reynolds number r v in the 
range 1.5 to 2.0 is necessary for convergence. Above a certain 
value of Reynolds number the iteration procedure would not converge 
for any value of r v . 
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Boundary-Layer Formulation 


The boundary-layer equations are expressed in terms of 
Gortler variables in order to reduce the variation of bound- 
ary-layer thickness in computational space. Thus the new co- 
ordinates s and n are introduced as 


s = x/c , n = y/6 


where x and y are distances measured along and perpendicular 
to the airfoil surface (or wake centerline), c is the airfoil 
chord length, and 6 is a thickness scale defined as 


6 = [ ( 2vc/u^ ) 

5 



u s ds] 


1/2 


The non-dimensional stream function F(s,n) is introduced, so that 
the x— component of velocity is given as 

u = u g F 1 (s ,n) 


where the prime denotes the n-derivative. The displacement 
thickness also is given in terms of as 


6 * = \ (l-u/u e )dy = 6(s)[n e - F(s ,n e )/F* (s,n e ) ] 


The boundary-layer momentum equation then becomes 

pt.t + ppt » + $ l - 0F» 2 = e 2 (pt || p»*) (5) 
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where 


3 


r2 du 

o s 

vc ds 


6 


2 


3 


1 


u s 6 2 /vc 

? u e 

(6 /vc) -p 

s 


du 

ds" 


e 


The displacement thickness is rescaled as 

A* = 6*/<5 - n e - P(s,n e )/P ' (s,n e ) (6) 

Let the subscript i denote the i— th grid point, as 
F i ( n ) = Fts^n). Then in these boundary-layer variables the 

interaction condition ( *0 becomes 

M 

Pi (n e ) - (r v Cil 5 i )A I = 1 * 3 l 1 C U 4 ! ‘ r v c ii 4 i (7> 


where 

The interacting boundary- layer problem is defined by Eqs. (5)-(7). 

In this formulation the Reynolds number R M occurs only through 

1N - 1/2 

the variables 6 and 6*, both proportional to R N 

It may be noted that u g = u s for the direct problem, in 
which case 3^ = 3; also in the similarity case F and F* are 
independent of s so that the right side of Eq. (5) vanishes, 
and the classical Falkner-Skan equation results. For the inter- 
action problem, u e (and hence 3 X ) is unknown and must be deter- 
mined along with the velocity profile F£(n). 
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Numerical Procedure 


The computer program is based on a second-order accurate 
finite-difference representation of Eqs. (5)- (7). The stream 
function was chosen as primary variable to facilitate treatment 
of as suggested by Eq . (6). Four-point centered differencing 

is used in the n-variable; typical terms are 

F " ,( ’Vi/2 ) - (P J + 1 - + 3Fj_ 1 - Fj_ 2 )/(in) 3 

p ” (n J-l/2 ) - <Vl - F j - Vi + Fj. 2 )/2(4n) J 
P ' (n J-l/2 ) =(F j - F j-1 )/4 ' 1 
P(n 3-l/2 ) * (F j + F j-1>/ 2 

On the other hand, backward-differencing is used in the s-variable. 
Either two— point or three— point differencing is available as an 
option (NX = 2 or NX = 3). For two-point differencing 

( H>i • (F i - F i-p/ 4s 

while for three-point differencing 
3 F 

= (3F i - + P 1 ^ 2 )/2As 

For the direct problem, on n=0 the no-slip conditions are 

F = F* =0, while on the wake centerline F = F" *» 0. Also since 

= u„ 
e s 

F* = 1 at n = n e (8) 

In this case the coefficients of the difference equations form 
a banded matrix having four non-vanishing diagonals, two of which 
lie below the principal diagonal. This structure is unchanged by 
Newton-Raphson iteration and the solution is obtained by a variant 


of the Thomas algorithm, described below. 

For the interaction problem, the edge condition (8) is 
replaced by condition (7), which can be written as 

FjJ + KA* = R (9) 

Eliminating A* by use of (6) gives 

< P N> 2 + K <1 N P N - V - RP N 

Since F” = 0 we can discretize F^ as (F^ - F ^_^)/h. Thus the 
(N+l) st equation is 


where 


A N+1 F N-1 + b m+i"m 0 


A N+1 ” P N-1 - P N - h(Kn N - R> 


Vi ■ - Kh - A N+1 


Linearizing gives the final form for Newton iteration: 


where 6F. 

j 

and 


A N+1 6F N-1 + B N+1 6F N “ R N+1 (10) 

is the change of F^ from one iteration to the next. 


A* « 

a N+1 

A N+1 “ f n 

+ 

F N-1 

B# = 

C N+1 

b n+i + f n 

- 

F N-1 

R# = 

n N+l 

“ A N+1 F N-1 

+ 

b n+i f n 


The interaction problem introduces an additional unknown, 
to permit use of the same programming as for the direct 
problem, the element Fj for j - N refers to the stream function 
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while the element F, J+ ^ refers to B-^. The matrix structure is 
still banded, but in addition its right-most column is non-zero. 
Thus, the complete system of equations for the Newton iterates 
has the structure given by (11) below. [Note that P = 0 for the 
direct problem, while represents the error in the original 
difference equations,] The subdiagonal elements are eliminated 
line-by-line, proceding from top to bottom; the system of re- 
duced equations is then solved from bottom to top. Convergence 
of the Newton iterates is rapid. For a boundary layer on a 
solid surface, convergence to | 6 F | max < 10 ^ is achieved in four 
or five iterations with the direct problem, and perhaps double 
that number for the interaction problem. On passing into the 
symmetric wake, the iteration count rose to about 20 or 25, 
falling slowly to the original four or five farther downstream. 



For the symmetric wake, the system of equations (11) actually 

begins with an equation for 6F 0 , the equation for 6F 0 also contains 
* d * j 

tne element B^, and that for 6F^ contains as well. For the solid 

wall (located at = 0, i = 1.5), the no-slip condition gives 

F^ “ F2 = 0 so that the first unknown is F^, as indicated in (11). 
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For the symmetric wake, F" vanishes instead of F', so that F^^ and 
F 2 are not known explicitly. Hence the differential equation (5) 
is solved on n = 0 as well as for n > 0. Using the four— point 
differences listed above we have 

1 

A 2 P 0 + b 2 f x + D 2 P 2 + s 2 P 3 + P 6 X - R 

f 

(Note P = 0 for the direct problem, while R = 0 for the inter- 
action problem). From the symmetry of the wake flow 


Hence the centerline equation reduces to the form 
(D 2 - B 2 )F 2 + (S 2 - A 2 )F 3 + P3 1 = R 
while for the Newton iterates 

D*6F 2 + S 2 6F 3 + ?$ 1 = R* 

The solution procedure is unchanged. 

The solution of the boundary-layer equations for a single 
streamwise station is carried out in a subroutine. The Hilbert 
integral (sum) in Eq. (7) is evaluated in the main program, along 
with mesh control and overall iteration of the boundary-layer sweeps. 
The mesh thickness in the n-coordinate is varied to account for the 
exaggerated growth of the boundary-layer thickness when separation 
occurs. A was used to control the mesh thickness, according to 
the equation 

(N - N Q ) An » A - A 0 

where the subscript 0 denotes values with which each boundary- 
layer sweep are initiated at the leading edge (stagnation point). 
This scheme permits the shear layer to be well separated from the 
surface and yet retains the resolution in n needed to allow for 
the strong shear at the dividing streamline. Mesh-thickness 
control was found to be essential for the success of the method. 
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Veldman also used a variable thickness grid, although he used a 
different technique. 

Some special features of the program are summarized below: 

O Both direct and interaction modes are allowed. The 
direct mode is used at the leading edge, a mixed mode 

is used for 0 < x < x 1 , and full interaction is used for 
x >_ x r 

OThe Rehyner/Flugge-L’dtz (Ref. 15) technique (FLARE) is used 
to stabilize reversed-f low computations. 

©Newton iteration is used for rapid convergence at each 
streamwise station. 

©A non-uniform grid is allowed in the streamwise direction, 
with optional subdivision to improve resolution when needed. 

©Optional two-point or three-point backward - differencing 
of the boundary-layer equations is available in the stream- 
wise coordinate. 

©Variable mesh thickness in the n - coordinate is included 
to account for boundary-layer growth when separation occurs. 
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DISCUSSION OF RESULTS 


Test Cases 

Two test cases that exhibit airfoil-like features were run as 
a check on the program: (1) the flat plate with wake at finite 

Reynolds number, and (2) the Carter-Wornom dent in an infinite plate. 
Both of these cases were computed by Veldman (Ref. 1), and the flat 
plate was computed by Werle and Verdon (Ref. 16), using a semi-in- 
verse iteration scheme due to Carter (Ref. 17). 

The flat-plate calculations were carried out for a Reynolds 
number of 10^ for both uniform and non-uniform grids (x-only), the 
latter corresponding to Carlson’s airfoil grid referred to above. 

For the uniform grid, the results compared very well with those of 
Veldman and of Werle and Verdon. With the non-uniform grid, it 
was necessary to refine the wake grid-point distribution to achieve 
satisfactory accuracy downstream of the trailing edge. This point 
will be emphasized later in the airfoil study. Convergence was 
rapid for this unseparated flow, with successive iterates for 6 * 
agreeing to five significant figures after ten iterations. 

Computations for the Carter-Wornom dent were made for R^ = 80,000, 
with good agreement with Veldman’ s results. In this case, variable 
grid thickness was necessary for accuracy owing to the extreme vari- 
ation of boundary-layer thickness at separation. Convergence was 
slower than for the unseparated flat-plate flow, requiring *13 iter- 
ations for convergence of 6# to five significant figures. 

These test cases incorporate singly two difficult features of 
viscous interaction: flow separation and wakes. These features 

combine in flow past airfoils, and it will be seen that convergence 
is much slower when flow reversal occurs in wakes. For comparison, 
in the results below for the NACA 0012 airfoil at a Reynolds number 
of 100,000, convergence to four significant figures required 6 0 
iterations. 
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Airfoil Results 


Interacting flow computations were carried out for two 
symmetric airfoils at zero incidence: the NACA 0012 Airfoil 

and the NACA 663-018 airfoil (see Figure 2). The 0012 section 
shape is characterized by adverse pressure gradient over most of 
the chord, whereas the 663 -OI 8 shape has minimum pressure at 60 
percent chord with adverse pressure gradient aft of that point. 
Hence these two airfoil sections represent a wide range of flow 
conditions. Computed results are presented in Figures 3-8 for 
the 0012 airfoil and in Figures 9-11 for the 663 -OI 8 airfoil. The 
effects of computational-grid parameters are indicated in Figures 
3-5. Of these parameters, the point of initiation of interaction 
and grid resolution are the two most important. 

For the non-interacting boundary layer on the 0012 airfoil, 
laminar separation occurs at about 60 percent chord; hence inter- 
action must begin ahead of that point. Various locations of the 
point (x^) of initiation of full interaction were attempted from 
= 0*19 to x 1 = 0.52. Converged results for these extremes are 
shown in Figure 3» in terms of the computed 6 # versus chordwise 
location from leading edge (x = 0 ) to about one chord length 
beyond the trailing edge (x/c = 1 ). The results for x.^ > 0.19 
are unsatisfactory since the interacting and non-interacting 
flow solutions deviate beyond that point. Consequently the re- 
maining results are for x 1 = 0.19. This value is satisfactory 
for both lower and higher values of Reynolds number, since the 
separation point is farther aft at lower Reynolds number, and 
the length of interaction is reduced at higher Reynolds numbers. 

The effect of grid resolution is shown in Figure 4. As 
mentioned earlier, Carlson's program was used to generate the 
non-interacting potential-flow solution, and the interacting 
flow solution is then computed on the same non-uniform grid. The 
results shown here correspond to an x-grid size of about three per- 
cent chord on the airfoil, but considerably larger in the wake. 



The actual wake-grid points for Carlson's grid are indicated 
by vertical lines in the figures. A fine grid was generated by 
interpolating points in the wake such that a nearly uniform 
grid (Ax - 0.03c) is obtained from leading edge to the wake cut- 
off position Xg, where the boundary-layer solution is stopped. 

Figure 4 shows that use of Carlson's grid leads to very large 
errors in the wake and, owing to upstream interaction, to signifi- 
cant errors on the airfoil as well. Consequently, only fine-grid 
results are presented in the subsequent figures. Computation time 
is increased greatly for the fine grid, since the number of points 
is doubled in each successive grid interval of the original coarse 
grid aft of the trailing edge. Thus for N points along the wake 
centerline in the original grid, 2(2 N -1) points are used in the 
fine grid. In particular, N = 5 for the cases shown in Figure 4 
(*2 ° 2.8c); thus 62 points are distributed along the wake center- 
line for the fine grid, versus the original five points. In both 
cases, 33 grid points were distributed along the airfoil. 

The effect of wake-cutoff position X 2 (where the boundary- 
layer solution terminates) is shown in Figure 5, corresponding 
to two wake lengths of about one and two chord lengths. The 
shorter wake appears to be sufficiently accurate, especially on 
the airfoil itself. Shorter values were attempted (X 2 = 1.4c), 
but were not satisfactory since the reversed-f low "bubble" 
does not close in that distance (for the NACA 0012 airfoil). 

Figure 6 compares results of computations using 2-point and 
3-point backward differencing for x-derivatives (MX = 2 and 3, 
respectively). The more accurate 3-point differencing yields some- 
what smaller values of 6* in the separated-f low region. (The 
separation point is at x/c = 0.69 for NX ® 3, 0.67 for NX = 2.) 
Recall that the approximation of Rehyner and Flttgge-Lotz (Ref. 15) 
is used, leading to some inaccuracy in the reversed-flow region. 
Since most of the computations were made with NX = 2, and Figure 6 
indicates that they are reasonably accurate, the following Reynolds- 
number comparisons are for NX = 2. 
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The effect of Reynolds number on the distribution of skin 
friction is shown in Figure 7. At the lowest Reynolds number, 
curve A, R^ = 1000, the flow is completely unseparated and the 
skin-friction has the highest values. As the Reynolds number 
increases, separation first occurs at the trailing edge, moving 
forward to about 70 percent chord at = 10,000, curve B. The 
skin friction then takes its lowest (most negative) value at the 
trailing edge. At a higher Reynolds number, curve C, R^ = 100,000, 
the separation point has moved forward to about 35 percent chord, 
and the skin friction reaches its most negative value at about 
^0 percent chord, slowly decaying in magnitude toward the trailing 
edge as the thickness of the reversed-f low region increases on the 
airfoil. At this Reynolds number, the separation "bubble” does not 
close in the wake even for x ? = 2.8. Attempts to carry out compu- 
tations at a higher Reynolds number (R fJ = 130,000) failed to con- 
verge. This lack of convergence for large values of R^ is discussed 
below. 

The effect of Reynolds number on the pressure distributions is 
shown in Figure 0 for the same cases. Viscous effects lower the 
suction peak at about 12 percent chord for all three cases, and 
likewise eliminate the stagnation pressure recovery at the trail- 
ing edge in all cases (thus producing pressure drag). The effect 
of flow separation is most pronounced for R„ = 100,000 (curve C), 
which shows a nearly flat pressure distribution over the whole 
separated region, at least beyond the point of minimum skin-friction. 
This result agrees with the classical Helmholz-Kirchhof f free- 
streamline theory, recently combined with triple-deck theory by 
F. T. Smith (Ref. 18) for inviscid separated flow behind bluff 
bodies . 

Turning now to the NACA 66^-018 airfoil section, we expect 
certain differences in the viscous characteristics to be ap- 
parent. First, the airfoil is designed to have a favorable pres- 
sure gradient over the forward 60 percent of the chord length, 
followed by a rather severe rising pressure. Hence, the laminar 
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separation point should be nearly fixed in position.* Second, 
the relatively thick section combined with the rearward location 
of separation should produce a rather thick separation "bubble" 
that extends rather far into the wake. 

A few computations have been carried out for the 66^-018 air- 
foil at Reynolds numbers of 10,000 and ^0,000. In this case, the 
inviscid surface velocity distribution was taken from the tables 
of Abbott and Von Doenhoff (Ref. 19) and supplemented by use of 
the Hilbert integral for values in the wake. The grid was es- 
sentially uniform in the x-direction, with Ax = 0.05c, but with 
finer intervals near the leading edge. Full interaction was 
initiated at 30 percent chord, since separation is delayed to 
about 60 percent chord on this airfoil. The wake was computed 
only to one chord length downstream of the trailing edge, al- 
though the results indicate that it should be continued farther- 

Figure 9 shows the resulting skin-friction distributions. 
Separation occurs at about 55 percent chord at both Reynolds 
numbers compared with the 60 percent location deduced from the 
inviscid pressure distribution. The most negative skin friction 
values occur at 60 percent chord, followed by a slow decay to- 
ward the trailing edge. This trend is like that for the 0012 
airfoil at R^ = 100,000, corresponding to the thicker separated 
layer. 

The corresponding pressure distributions are shown in 
Figure 10. Their behavior is similar to the R^ = 100,000 case 
for the 0012 airfoil, with deviation from the Inviscid distri- 
bution occurring near the leading edge and with a suction mini- 
mum in the separated region. Unlike the 0012 case, the pressure 
tends to fall (suction rises) as the trailing edge is approached. 
However, the computed reversed-f low "bubbles" for both cases 


Note that at high Reynolds number, the boundary layer is turbu r 
lent and does not separate on the 66^-018 airfoil at zero 
incidence . 
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A and B in Figure 10 do not close up to the point at which 
the wake computations were terminated. Consequently, these 
results must be viewed with some reservations. 

Computations for the 66 3 ~0l8 airfoil at Reynolds number of 
100,000 were not convergent, corresponding to a similar failure 
for the 0012 airfoil at a higher Reynolds number. It is thought 
that this non-convergence may be caused by thinning of the 
separated shear layer as Reynolds number increases. This hy- 
pothesis is illustrated by Figure 11, which presents the boundary- 
layer velocity profile at the trailing edge for R^ = *10,000, in 
terms of the boundary-layer coordinate n = y/6 (see above Eq 5). 

At this Reynolds number, the nearly stagnant reversed-f low layer 
is about as thick as the shear layer above it. This comparison 
would be more pronounced at a higher Reynolds number. A velocity 
profile of this type is highly unstable, so that turbulence would 
develop under such conditions. Flow visualization studies of the 
6^ 3 — 018 airfoil by Mueller and Batill (Ref. 20) show laminar flow 
to the trailing edge at R^ = *10,000, with oscillatory vortical 
flow immediately downstream, whereas at Rj. = 130,000 the entire 
separated shear layer appeared to be turbulent. At R^ = *100,000 
no separation occurred (for zero incidence), probably because the 
boundary layer became turbulent upstream of the laminar separ- 
ation point. These comparisons suggest that the breakdown of the 
present laminar interaction computations coincides with conditions 
for which laminar flow ceases to exist in the separated shear layer. 
The conclusion then is that a reliable boundary layer-transition 
model is needed for extending the computations to moderately 
higher Reynolds number. 


CONCLUSIONS 


The computational results show that the present method for 
solving the laminar interacting boundary-layer equations is most 
successful at relatively low Reynolds number where the separated- 
flow region is small. Convergence is progressively more difficult 
as the Reynolds number is increased, and the method fails to con- 
verge at some value of R { j that depends on the airfoil shape. This 
breakdown of convergence appears to correspond to conditions for 
which the viscous shear layer is thin and well-separated from the 
surface. Such shear layers are highly unstable and transition to 
turbulence would be expected. Hence, it is concluded that the 
method works best where laminar flow occurs naturally. 

Although the present work was restricted to symmetrical 
flows, airfoils at angle-of attack could be treated by solving the 
boundary-layer equations (and interaction condition) separately on 
both sides of the airfoil. Similar computations have been carried 
out using the triple-deck model (Daniels, Ref. 21, Chow and Melnik 
Ref. 12, Mansfield, Ref. 22). However, the separated region would 
become thicker, and the breakdown of convergence would be expected 
to occur at lower Reynolds number. The remarks above concerning 
transition to turbulence apply here as well. The method easily 
extends to turbulent flow by incorporating a turbulence model in 
the bound ary -layer equations as already carried out by Burggraf 
(Ref. 23 ). However, at intermediate Reynolds-numbers where trans- 
ition occurs in the separated shear layer, a transition model 
would be required. While such models exist, they are not very 
reliable. Another area of interest is that of leading-edge sep- 
aration, for which all of the above remarks apply. In addition, 
the present method would need to be modified to permit full inter- 
action at (or very near to) the leading edge. In that case, the 
use of the Hilbert integral as interaction condition would no 
longer be appropriate; direct coupling of the boundary-layer com- 
putation with a full potential-flow solver would be recommended. 
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APPENDIX A 

EVALUATION OF HILBERT INTEGRAL ON IRREGULAR GRID 


Two discretized formulations of the Hilbert integral have 
been given by Veldman (Ref. 1*0, both for uniform grids. For 
our purposes, it is convenient to allow a non-uniform grid, 
and so one of his algorithms is generalized here for that 
situation. Denote the arbitrarily spaced coordinate values 
by Xj , j = 0,1,..., N+l, where Xq < < . . . < x^ + ^. Then the 

Hilbert integral in Eq. (1) can be replaced by a finite sum 
with each term centered in the grid intervals, as 


r x n+l 

.. d6_ dx ^ J 
dx x^-x *** i 


N 

E 

j=0 


rd6 -| 

ax J+l/2 


( x j +l” x j)/( x i”” x j + l/2 


) 


where x J+1/2 = (Xj+x j+ i)/ 2 , and we drop the superscript (*) for 
convenience. The limits on the integral Xq and x^ + ^ are assumed 
to be sufficiently wide that dS/dx is negligible beyond those 


points. Using central differencing 


r d6n 

L dx J j+l/2 


^ 6 j+l“ 6 J ^ X J + l“ X J ^ 


Hence 


J i = J ^ 0 (6 j + l" <S j )/(x i" X J+l/2 ) 


The sum in 
obtain 



can be reindexed and 


combined with that in 


N 

= - E (Aj/ 2 ) (x j + 1 -Xj_ 1 )/(x 1 -Xj__ 1/2 ) (x 1 -x J + 1/2 ) 

J 

- A 0 /(x r x 1/2 ) + A N+1 /( Xi - x N+1/2 ) 


to 
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In terms of the notation of Eq. (3.3), we have 
c iQ = -l/( Xi -x 1/2 ) 

c i,N+l = 1/(x i” x N+l/2 ) 

This formulation of the Hilbert integral is second-order 
accurate, consistent with the differencing of the boundary-layer 
equations. Note that the Cauchy principal value of the inte- 
gral is evaluated accurately by centering the integrand between 
adjacent grid points. For a uniform grid, this algorithm reduces 
to Veldman*s simpler version. 


28 



APPENDIX B 


MULTIGRID MODIFICATIONS TO A POTENTIAL-FLOW CODE 


The Multi-grid Concept * 

Suppose we want to solve the following finite-difference analog 
to some differential equation: 

^■h^h c F h (B-l) 

(The subscript h signifies a finite-difference representation 
on a grid with spacing h.) Further, let us assume that we have 
some approximation <f>^ to and let the difference between the 

two be i.e. 

*h " $h “ *h (B-2) 

Then, if is a linear operator (as it is in our case) 

^h^h = ^h $ h + ‘Oh = F h (B_3) 


or 

®fh ? h = F h “^h $ h = R h 


(B-lt) 


The idea of multigrid is to make use of the fact that we can 
make useful (and inexpensive) approximations to the above equa- 
tion (B-4) on a coarser grid, say of mesh size 2h. Thus using 
some relaxation technique, we solve the equation 


X * 


- T 2h R 

2h 2h i h K h 


(B-5 ) 


where 1^ is an operator which transfers to each point on the 
coarse grid the corresponding value of the residual R h on the 
fine grid. The solution is obtained more rapidly on the coarse 
grid than on the fine grid. Once we have obtained the solution 
of (B-5) (or at least a good approximation to it) on the coarse 


Reference 24 is an excellent treatment of this subject for 
general applications. 
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grid, we interpolate the correction $ 2h to the finer grid: 

*h = I 2h ? 2h 
h 

where I 2h is an Interpolation operator. Depending upon a 
number of factors (grid geometry, interpolation order), the in- 
terpolated value of $ h may or may not satisfy (B-4) closely 
enough to be used as a correction to If not, we simply 

relax (B-lJ) a few times and then use 


^h 

new 




(B-6) 


to get a better approximation to <J> h , which is the solution of 
the finite-difference equation of interest. In general, the multi- 
grid procedure uses a sequence of grids, rather than only 2 grids. 
To achieve fast convergence, the multigrid method relies on the 
fact that standard relaxation techniques converge rapidly during 
the initial iterations, but exhibit progressively slower conver- 
gence at later stages, especially for fine grids. This is be- 
cause relaxation techniques are efficient at removing high fre- 
quency-error components (i.e., errors with wavelengths on the 
order of the grid spacing h), but inefficient at removing low 
frequency-error components (i.e.» errors with wavelengths on the 
order of the overall grid size). The multigrid method circumvents 
this difficulty by using a series of coarser grids to eliminate 
the low frequency-error components. In effect, each of the coarse 
grids acts as a filter which removes error components with cor- 
responding wave lengths. Peaceman and Rachford (Ref. 25) 
formulated a method for solving Laplace's equation exactly in a 
finite number of steps, with one spectral component of the error 
removed in each step. For NxN grids their method requires of 

■p 

order N operations to obtain the exact solution, whereas the 

p 

multigrid method requires only of order N log N operations to 
reduce the error to a specified tolerance. 
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Application to Potential Flow Pa3t Airfoils 

Carlson's (Ref. 1) code has been modified to incorporate 
the multigrid solution procedure. For efficiency, the transonic 
code was simplified to allow only the incompressible case. The 
resulting code then solves the difference equation 


7 h* h ■ 0 

subject to 

( 3x> “ " [(slna + <'hy )/(oos ° + 

s s 

where the subscript s denotes the airfoil surface, and the 
subscripts x and y denote differentiation with respect to x 
and y > and 

4> h -*■ r h (0-a)/2ir as x 2 + y 2 -*- 00 

where tan 0 ** y/x and I* h ■ <j> h (x,0 + ) - <J>^(x,0"’)for x > x TE . 

Substitution of 4^ = $h + ^h al30Ve equations leads 

to the following residual equations for 

v h* h - "h 11 (B - 7a> 

subject to: 

( S ) s (? h v ) 8 " (? h } s “ R h 2) CB-7b) 

x y 

$ h - - f h (e-a)/2ir as x 2 + y 2 ^ 00 (B-7c) 

where f h = $ h (x,0+) - $ h (x,0“) for x > x TE . 


Here the residuals R^j^and R^^ 
residuals, respectively; 

R< “- - 7 

R^jp » sin a + ($ h ) s -[cos a + 


denote interior and boundary 


^h (B-8a) 

( ^h_ ) 8 3( S ) 8 (B-8b) 
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Coordinate Stretching 


The coordinate stretching scheme of Carlson has been pre- 
served in the multigrid routine. Carlson (Ref. 1) shrinks both 
the x and y directions onto the finite computational plane 
£ ,n according to 

x = + A 2 tanO(S + 5^/2] + A^anOU + Cj, ) 3 /2] forlx^Xj, 

x = £(a + b£ 2 ) for |x|<x^ 


and y = A x tan(7rn/2) 

where A^, A 2 , A^, a, b, x 

f ” 5x 8111(1 S “ 


and 


dn 

3y 


are constants. 


Then, denoting 


equations (B-7a) and (B-7b) can be written as 
f(f? 5 ) 5 + g(g$ n ) n = R (1) 


(B-9a) 



< f Vs 


<e*n>s 


r ( 2 ) 


(B-9b) 


where the subscript h on ? and R is understood. In the new 
coordinate system, the circulation correction ? is determined 
from: 


? a $(5,0 + ) 
The residual functions 


- ?(C,0~) for £>S te 
R^ 1- ^ and R^ 2 ^ become 


(B-10) 


R (1) - - g(g$ n ) n (B-lla) 

R (2) = sin a + (g? n ) s - [cos a +(« 5 > s ] (f£) s (B _ llb) 
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Finite Differences - Governing Equations 

Employing second order-accurate differences In £ and n* 
equation (B-9a) can be written as 


tg J g J-l/2 /(An) ] *ij-r C(f i+l/2 + f i-i/2 )f i /a, ( A ° + (g J+l/2 +g J-l/2)* 

. gj /(An) 2 ]^ + Cg j 6 j . 1 / 2 /(An ) 2 ]^ >J+1 = 

-Cf ± /(A5) ^ f i+1 / 2 ?i+lj " (f i+l/2 + f i-l/2 )(1 ” 1/a, ^lj 


where as In Carlson (Ref. 1), to Is a relaxation factor and 
+ denotes new values. The above represents a tridiagonal system 
of equations for the values of (fi^, on a given column, 1. 

Equation (B-12) is solved on columns, sweeping from upstream to down 
stream. 

The residual function must also be replaced by its 

finite-difference form. Using second order-accurate differences 
in 5 and n, equation (B-lla) is represented as: 

R iJ ^ “ “ Cf i /(A5)2][ - f i+l/2^i+l,j“ (f i+l/2 + f i-l/2^i,J 

f i-l/2^i-l,J] - Cgj/(An) 2 ][gj + 1/2 ^ij + i " g j-l/2 * 


0 



$ i,J-l } “ g J+l/2 $ i,j ] 


(B-13) 


Finite Differences - Boundary Conditions 

As in Carlson's original program, the residual tangency 
condition, equation (B-lla) is used to create fictitious values 
of <j> inside the airfoil, allowing equation (B-12) to be used as 
it stands at points immediately outside the airfoil. We expand 
(4> n ) g 501(1 (^) s in Taylor series about fictitious points in- 

side the airfoil and replace the resulting derivatives of <J> by 
finite differences yielding second-order accurate results. 
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for the upper surface 

* 2(^) 2 /e s [3An - 2 (n s - n JA _ a )r 

’{ R iu ) + E s [(1,? i,JA ' ? l,jA-H )/2in 
' (n s • n jA-l H2? l,jA - ? i,jA + l )/(4n)2] 

( to ) s r 5 ^1+1,JA-1 ' ^i-l,JA-l* /2AS * (n s ' n JA-l'’ 

* (? 1+1,JA " ? l+l,jA-l * ? 1-1,JA + , l-l,jA-l )/2A5an:1 } 
and for the lower surface: 

^i,JA+l = -2(An) /g s C3An + 2(n s - a+1> ] R ±? ^ 

+ g s + $i,jA-l )/2An ~ ( n s ~ ^JA+l^^iJA 

‘ $ i,JA-l )/(An)2] ’ ( &s f s C($ i+l,JA+l “ $ i-lJA+l )/2Ae 

+(n s” n jA+l ) ($ i+l,jA+l“ $ l+l,JA“ $ l-l,jA+l +$ i-lJA )/2ACAn | (B “ 15) 

where the Index jA denotes the closest point outside the air- 
foil boundary , and R^ u and R^ denote residuals for upper and 
lower surfaces. 

Equation (B-14) or (B-15) is first solved on column i 
using old values from the previous relaxation sweep. Then 
equations (B-12) are solved for new values exterior to the air- 
foil along the column i. The new values of $ so obtained are 
used in equations ( B— 14 ) or (B-15) to generate new values for 
the fictitious interior points. Thus the interior points are 
based on new values at station i and old values at other 
stations as are all points external to the airfoil. 

The values of r£ 2) and R^ 2) are given in finite difference 
form by 
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R iu } " Sin a -<S>8 C0s a + + ^i,JA - *i,JA + l )/2AT1 

+ (n s ~ n JA-l )( ^i,jA-l “ 2 ^i,JA + $ iJA+l )/(An)2 ^ 
~*ax*s f s^i+l,JA-l “ $ i-l,JA-*l )/2A? 

+ (n s ** n JA-l^i+l,JA~*i+l,JA-l “*i-l,JA + ^i-1, j A-l )/2A5An3 

R lV m Sin a “ l$x>» 008 a + 6i (3 ^i„JA+l " ^i,JA + *i,JA-l )/2An 
+ (n s “ n JA+l )( ^i,JA+l “ 2$ l,jA + ^i,jA-l )/(An) J 

( dx^ f s C ^i+l,JA+l “ $ i-l,JA+l )/2A? 

* (n s " n JA+l )( ^i+l,JA+l “ *i+l,JA “ ^i-lJA+1 + ^i-l,JA )/2A?An3 

Computational Procedure and Results 

1. Three grid spaclngs are used: fine, medium and coarse. 

2. An initial approximation on the fine grid is obtained by 
first obtaining an approximation by relaxing on the coarse 
grid and interpolating these results to the medium grid, 
where they serve as an initial approximation. Relaxation 
is again performed on this grid, and these results are 
interpolated to the fine grid to be used as its starting 
approx imat ion. 

3. The residuals are calculated on the fine grid using (B-8) 
and are transferred to the coarse grids. 

4. A few relaxation sweeps (three are used here) are used on the 
coarse grid to solve (B-7) (approximately) on that grid. 

5. Using quadratic interpolation the coarse grid values of $ 
are transferred to the medium grid to serve as an initial ap- 
proximation to (B-7) there. At the same time, the residuals 
are transferred to the medium grid. 
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6. (B-7) is relaxed a few times on the medium grid and the 
results Interpolated to the fine grid. Three relaxations 
suffice here . 

7. In the present problem, it was found necessary to relax 
three times on the finest grid, after which the new approxi- 
mation to <J> was found using (B-6). 

8. On the fine grid, we relax cf> three times and check the 
convergence tolerance: max|A<J>| < e. 

If convergence has not been achieved, new residuals are 
calculated and Steps 3 to 8 are repeated until convergence to 
the desired tolerance is achieved. Optimization of the multigrid 
cycle was not attempted. 

Multigrid computations were made at both zero incidence and 
at an incidence angle of five degrees. In the former case, the 
method gave the expected improvement in efficiency, but for the 
latter case there was no significant improvement over the 
original code. The Interior flow computations appeared to be- 
have well, but the circulation about the airfoil converged very 
slowly. It is possible that an error exists in the program 
relative to the treatment of circulation on the various subgrids, 
though a considerable effort to find one was not successful. 

To illustrate the method, we take the case of an NACA 0012 
airfoil at zero angle of attack . For comparison we list the CPU 
time on an Amdahl 470 V/6-II computer and the equivalent number 
of fine-grid iterations used (i.e., iteration on the coarse grid 
of spacing 2h is 1/4 of an iteration on the fine grid) . The same 
value of relaxation factor was used throughout. The convergence 
tolerance was set at 10 . 


Finest Grid 

No. of 
Grids 

Original Code 

Multigrid 

Code 

CPU(sec) 

Work Units 

CPU (sec) 

Work Units 

25 x 13 

2 


20 

2.2 

46 

49 x 25 

3 

19 -5 

187 

5 .8 

cn 

in 

j 97 x 49 

1 

3 

59 -0 

114 

30 .6 

' 

1 

71 
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Figure 1. Summary of Triple-Deck Results for Incompressible 
Trailing-Edge Flow for a Flat Plate at Zero Incidence. 
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Figure 2. Configuration of Two NACA Airfoils. 
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Figure 3. Effect of Point of Initiation of Interaction on 
Displacement Thickness, for NACA 0012 Airfoil. 

R n = 10,000, NX = 2, Carlson Grid, X ? = 1.9. 

(11) X x = 0.52; (5) X ]L = 0.19. 
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Figure Effect of Grid Distribution on Displacement Thickness 
for NACA 0012 Airfoil. R N = 10,000; NX = 2; * 0.19j X 2 = 

(9) Carlson Grid; (1*0 Fine Grid. 






YSC 


Figure 5. Effect of Downstream Extent of Grid on Displacement 
Thickness for NACA 0012 Airfoil. - 10,000; NX = 2; 
Fine-Grid Wake; » 0.19. 

UNO X 2 » 2.8; (18) X 2 » 1.9. 
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Figure 6. Displacement-Thickness Results for Two-point and 
Three-Point Integration Schemes for NACA 0012 Airfoil. 

R n » 10,000; X 1 ■ 0.19, X 2 ■ 1.9. 

(16) NX = 3; (18) NX - 2. 






Figure 7. Effect of Reynolds Number on Distribution of Skin- 
Friction Coefficient for NACA 0012 Airfoil. 

NX = 2; Fine-Grid Wake; X 1 - 0.19, X 2 - 2.8 
(A)B n = 1000; (B) R n » 10,000; (c) R N = 100,000. 
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Figure 8. Effect of Reynolds Number on Distribution of Pressure 
Coefficient for NACA 0012 Airfoil. 

NX = 2; Fine-Grid Wake; X 1 = 0.19; X 2 = 2.8. 

(A) R n = 1000, (B) R n = 10,000; (C) R N = 100,000. 
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Figure 9. Effect of Reynolds Number on Distribution 

of Skin-Friction Coefficient for NACA 66^-018 Airfoil 
NX = 2; Fine-Grid Wake; X x - 0.30; X 2 - 2.0. 

(A) R m = 10,000; (B) R„ * 40,000. 






Figure 10. Effect of Reynolds Number on Distribution of 
Pressure Coefficient for NACA 66^-018 Airfoil. NX «* 2; 
Fine-Grid Wake; = 0.30; X 2 = 2.0. 

(A) R^j = 10,000; (B) ^ = 40,000. 
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Figure 11. Velocity Profile in Boundary Layer 
at Trailing Edge of NACA 66^-018 Airfoil at 
R n ® 40,000. 
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